This repository has been archived by the owner on Oct 14, 2023. It is now read-only.
/
Exploring the New Horizons launch.mystnb
169 lines (122 loc) · 4.06 KB
/
Exploring the New Horizons launch.mystnb
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
---
jupytext:
text_representation:
extension: .mystnb
format_name: myst
format_version: 0.13
jupytext_version: 1.14.0
kernelspec:
display_name: Python 3
language: python
name: python3
---
# New Horizons launch and trajectory
Main data source: Guo & Farquhar "New Horizons Mission Design" http://www.boulder.swri.edu/pkb/ssr/ssr-mission-design.pdf
```{code-cell}
from astropy import time
from astropy import units as u
from poliastro import iod
from poliastro.bodies import Sun, Earth, Jupiter
from poliastro.ephem import Ephem
from poliastro.frames import Planes
from poliastro.plotting import StaticOrbitPlotter
from poliastro.twobody import Orbit
from poliastro.util import norm
```
## Parking orbit
Quoting from "New Horizons Mission Design":
> It was first inserted into an elliptical Earth parking orbit
of **perigee altitude 165 km** and **apogee altitude 215 km**. [Emphasis mine]
```{code-cell}
r_p = Earth.R + 165 * u.km
r_a = Earth.R + 215 * u.km
a_parking = (r_p + r_a) / 2
ecc_parking = 1 - r_p / a_parking
parking = Orbit.from_classical(
Earth,
a_parking,
ecc_parking,
0 * u.deg,
0 * u.deg,
0 * u.deg,
0 * u.deg, # We don't mind
time.Time("2006-01-19", scale="utc"),
)
print(parking.v)
parking.plot()
```
## Hyperbolic exit
Hyperbolic excess velocity:
$$ v_{\infty}^2 = \frac{\mu}{-a} = 2 \varepsilon = C_3 $$
Relation between orbital velocity $v$, local escape velocity $v_e$ and hyperbolic excess velocity $v_{\infty}$:
$$ v^2 = v_e^2 + v_{\infty}^2 $$
### Option a): Insert $C_3$ from report, check $v_e$ at parking perigee:
```{code-cell}
C_3_A = 157.6561 * u.km**2 / u.s**2 # Designed
a_exit = -(Earth.k / C_3_A).to(u.km)
ecc_exit = 1 - r_p / a_exit
exit = Orbit.from_classical(
Earth,
a_exit,
ecc_exit,
0 * u.deg,
0 * u.deg,
0 * u.deg,
0 * u.deg, # We don't mind
time.Time("2006-01-19", scale="utc"),
)
norm(exit.v).to(u.km / u.s)
```
Quoting "New Horizons Mission Design":
> After a short coast in the parking orbit, the spacecraft was then injected into
the desired heliocentric orbit by the Centaur second stage and Star 48B third
stage. At the Star 48B burnout, the New Horizons spacecraft reached the highest
Earth departure speed, **estimated at 16.2 km/s**, becoming the fastest spacecraft
ever launched from Earth. [Emphasis mine]
```{code-cell}
v_estimated = 16.2 * u.km / u.s
print(
"Relative error of {:.2f} %".format(
(norm(exit.v) - v_estimated) / v_estimated * 100
)
)
```
So it stays within the same order of magnitude. Which is reasonable, because real life burns are not instantaneous.
```{code-cell}
:tags: [nbsphinx-thumbnail]
from matplotlib import pyplot as plt
fig, ax = plt.subplots(figsize=(8, 8))
op = StaticOrbitPlotter(ax=ax)
op.plot(parking)
op.plot(exit)
ax.set_xlim(-8000, 8000)
ax.set_ylim(-20000, 20000)
```
### Option b): Compute $v_{\infty}$ using the Jupyter flyby
According to Wikipedia, the closest approach occurred at 05:43:40 UTC. We can use this data to compute the solution of the Lambert problem between the Earth and Jupiter:
```{code-cell}
nh_date = time.Time("2006-01-19 19:00", scale="utc").tdb
nh_flyby_date = time.Time("2007-02-28 05:43:40", scale="utc").tdb
nh_tof = nh_flyby_date - nh_date
nh_r_0, v_earth = Ephem.from_body(Earth, nh_date).rv(nh_date)
nh_r_f, v_jup = Ephem.from_body(Jupiter, nh_flyby_date).rv(nh_flyby_date)
nh_v_0, nh_v_f = iod.lambert(Sun.k, nh_r_0, nh_r_f, nh_tof)
```
The hyperbolic excess velocity is measured with respect to the Earth:
```{code-cell}
C_3_lambert = (norm(nh_v_0 - v_earth)).to(u.km / u.s) ** 2
C_3_lambert
```
```{code-cell}
print("Relative error of {:.2f} %".format((C_3_lambert - C_3_A) / C_3_A * 100))
```
Which again, stays within the same order of magnitude of the figure given to the Guo & Farquhar report.
+++
## From Earth to Jupiter
```{code-cell}
nh = Orbit.from_vectors(Sun, nh_r_0, nh_v_0, nh_date)
op = StaticOrbitPlotter(plane=Planes.EARTH_ECLIPTIC)
op.plot_body_orbit(Jupiter, nh_flyby_date)
op.plot_body_orbit(Earth, nh_date)
op.plot(nh, label="New Horizons", color="k")
```